function partial_x = Partial_derivative_x(u, boundary_condition, differencing_scheme, h, sigma)
% partial_x = Partial_derivative_x(u, boundary_condition, differencing_scheme, h, sigma)
% 
% Partial_derivative_x calculates the partial derivative of a function 
% u : R^2 -> R (a 2D matrix) with respect to x (the first coordinate).
%
% input:    - u : 2D image with scaler values
%           - boundary_condition    (DEFAULT = 'constant')
%           - differencing_scheme   (DEFAULT = 'central')
%           - h                     (DEFAULT = 1)
%           - sigma                 (DEFAULT = 0)
%
% output:   - partial_x : a 2D matrix with the values of the partial 
%             derivative
%
% See also: Gradient
%
% 2013-07-25 --------------------------------------------------------------
% 
% Matthias J. Ehrhardt
% CMIC, University College London, UK
% 
% matthias.ehrhardt.11@ucl.ac.uk
% 
% -------------------------------------------------------------------------
    
    if nargin < 2, boundary_condition = 'constant'; warning('boundary_condition not chosen. DEFAULT = %s', boundary_condition); end;    
    if nargin < 3, differencing_scheme = 'central'; warning('differencing_scheme not chosen. DEFAULT = %s', differencing_scheme); end;
    if nargin < 4, h = 1; warning('h for discretization not chosen. DEFAULT = %3.2e', h); end;
    if nargin < 5, sigma = 0; warning('sigma for smoothing not chosen. DEFAULT = %3.2e', sigma); end;  
    
    switch differencing_scheme
        case {'forward'} % Forward differencing scheme
            c_derivative = 1./h * [0 -1 1]';
        case {'backward'} % Backward differencing scheme
            c_derivative = 1./h * [-1 1 0]';
        case {'central', '3Point'} % Central differencing scheme
            c_derivative = 1./(2*h) * [-1 0 1]';         
        otherwise
            error('Choice of differencing scheme unavailable')
    end

    [I, J] = size(u);

    if sigma ~= 0
        c_smoothing = fspecial('Gaussian', [max(3, round(2 * sigma)), 1], sigma);
        c_derivative = imfilter(c_derivative, c_smoothing', 'full');

        boundary_length_smoothing = round(0.5 * (length(c_smoothing)-1));
        J = (1 : J) + boundary_length_smoothing;
    end
    
    boundary_length_derivative = round(0.5 * (length(c_derivative)-1));
    
    I = (1 : I) + boundary_length_derivative;
    
    U = Extent_x(u, boundary_condition, boundary_length_derivative);
    U = imfilter(U, c_derivative);

    partial_x = U(I,:);
    
    if sigma ~= 0
    	U = Extent_x(partial_x', boundary_condition, boundary_length_smoothing)';
        U = imfilter(U, c_smoothing);

        partial_x = U(:,J);
    end   
end